% Simulate price data
T = 120; % Number of periods
regime = [ones(1, 50), 2 * ones(1, 50), ones(1, 50), 2 * ones(1, 50)]; % Simulated regimes
price_changes = zeros(T, 1);
reaction_window = 3; % Analyze 5 periods after jump
price_reactions = NaN(length(jump), reaction_window);
% Generate price changes based on regime
%for t = 2:T
%    if regime(t) == 1
%        price_changes(t) = 0.02 + randn() * 0.1; % Normal regime
%    else
%        price_changes(t) = 0.05 + randn() * 0.5; % High-volatility regime
%    end
%end

% Cumulative price levels
%prices = cumsum(price_changes);
p = cumsum(P);
prices=p';

% Plot simulated data
figure;
plot(prices);
title('Simulated Prices with Regimes');
xlabel('Time');
ylabel('Price');
grid on;

% Define Markov Switching Model
num_states = 2; % Two regimes (normal, high-volatility)
msm = msVAR(price_changes, num_states);

% Estimate the model parameters
[est_params, states] = msm.estimate();

% Display estimated parameters
disp('Estimated Parameters:');
disp(est_params);

% Plot inferred states
figure;
plot(states, 'o-');
title('Inferred Regimes from Markov Switching Model');
xlabel('Time');
ylabel('State');
grid on;

% Identify jumps
jump_threshold = 0.5; % Threshold for defining a jump
jumps = abs(price_changes) > jump_threshold;

% Analyze inferred regimes after jumps
post_jump_states = NaN(sum(jumps), 5); % Matrix to store states post-jump

jump_indices = find(jumps);
for i = 1:length(jump_indices)
    idx = jump_indices(i);
    if idx + 5 <= T
        post_jump_states(i, :) = states(idx:idx+4); % Store 5 periods after jump
    end
end

% Calculate regime frequencies after jumps
regime_frequencies = sum(post_jump_states == 2, 1) / size(post_jump_states, 1);

% Display frequencies
disp('Proportion of High-Volatility Regime (Post-Jump):');
disp(regime_frequencies);

% Plot regime frequencies
figure;
bar(1:5, regime_frequencies);
title('Proportion of High-Volatility Regime After Jumps');
xlabel('Periods After Jump');
ylabel('Frequency');
grid on;

